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Abstract 

We  show  how  a  Bayesian  analysis  of  a  fertility  model  incorporating  many  of  the 
previously  suggested  models  can  account  for  uncertainty  about  which  fertility  model 
provides  the  best  approximation  in  any  given  trial.  We  also  show  how  uncertainty 
about  anomalies  such  as  outliers  and  fertility  jumps  can  be  accounted  for.  We  argue 
that  this  is  preferable  to  conditioning  on  an  “appropriate”  model,  and  show  by  examples 
how  accounting  for  such  possible  anomalies  can  both  influence  support  for  a  particular 
fertility  model  and  reduce  the  dependence  of  treatment  estimates  on  the  choice  of 
fertility  model. 
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1  Introduction 

The  object  of  most  agricultural  field  trials  is  to  assess  the  effect  of  treatment  on  yields.  This 
must  be  done  in  the  presence  of  unknown  fertility  trends  through  the  field.  Many  stochastic 
models  for  these  fertility  trends  have  been  proposed  and  shown  to  be  efficient  competitors 
to  classical  blocking  designs  (see  Section  2).  In  this  paper  we  shall  restrict  attention  to 
trials  where  plot  fertilities  are  correlated  in  only  one  direction.  In  the  case  of  a  rectangular 
lattice  of  plots  this  is  often  assumed  when  the  long,  thin  nature  of  the  plots  makes  those 
with  common  long  edges  spatially  close  compared  to  those  with  common  short  edges. 

When  using  a  stochastic  model  for  field  fertility  in  the  analysis  of  a  particular  trial  two 
questions  must  be  addressed.  The  first  is  the  choice  of  fertility  model,  and  the  second  is  the 
treatment  of  anomalies  such  as  outliers  and  fertility  jumps,  i.e.  sudden  level  shifts  in  fertility. 
Even  when  several  such  possibilities  are  entertained,  it  is  not  uncommon  to  condition  on  a 
single  model  and  a  single  set  of  anomalies,  which  are  determined  by  model  diagnostics;  sec 
for  e.xample  Cullis,  McGilchrist  and  Gleeson  (1989)  and  Martin  (1990).  Here  we  propose 
to  account  for  uncertainty  about  the  fertility  model  and  possible  anomalies  rather  than  to 
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condition  on  an  “appropriate”  model.  This  will  be  particularly  important  when  several 
choices  are  almost  equally  appropriate  (see  Section  4.3). 

We  first  define  a  general  but  simple  state-space  model  for  agricultural  field  trials,  which 
turns  out  to  include  many  of  the  models  in  current  use  (Section  2).  We  then  show  how 
this  framework  allows  us  to  model  outliers  and  fertility  jumps  quite  easily  (Section  3).  The 
importance  of  doing  this  and  of  accounting  fully  for  the  associated  uncertainty  is  illustrated 
in  several  examples  (Section  4). 

2  A  Fertility  Model 

2.1  Model  Definition 

We  suppose  that  the  trial  consists  of  one  or  more  blocks  of  plots.  Here  a  block  refers  to  a 
row  of  adjacent  plots,  with  plots  in  different  blocks  cissumed  to  bear  little  relationship  to 
each  other.  The  plots  are  numbered  1  through  noplots  along  the  row  within  each  block,  and 
block  by  block.  Plot  t  is  therefore  at  the  beginning  of  a  new  block  if  /  =  1  or  plots  t  and 
t  —  1  are  in  different  blocks. 

We  assume  that  the  data  can  be  decomposed  additively; 

Y  =  Dt-\-F-\-c  (1) 

where  Y  is  the  noplots-vector  of  yields,  t  is  the  notrcats- vector  of  treatment  effects  (which 
may  include  other  covariates),  D  is  the  corresponding  design  matrix,  F  is  the  noplots-vector 
of  fertility  levels,  and  e  is  the  noplots-vector  of  measurement  errors. 

We  model  F  recursively  to  correspond  to  walking  through  the  plot^  from  plot  1  to  plot 
noplots,  with  t  therefore  indexing  time  as  well  as  the  plots.  Starting  from  the  assumption 
that  the  fertility  is  approximately  locally  linear,  we  suppose  that  the  fertility  level  at  plot 
t  not  only  depends  on  the  fertility  level  at  the  previous  plot,  but  also  on  some  measure  of 
the  rate  of  increase  in  the  fertility  at  this  plot.  We  define  G'(_i  to  be  this  fertility  gradient 
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between  the  fertility  level  at  plot  t,  namely  Ff,  and  that  at  plot  i  —  \  ,  namely  F(_i. 

More  precisely,  we  shall  model  the  plot  fertilities  conditional  on  the  parameters  Ai,A2, 
and  where  0  <  A2  <  Ai  <  1  and  cr^ad  —  0  recursively  as  follows, 

case  1  :  same  block; 


(S),' 

case  2  :  new  block; 
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where  we  use  the  notation  (f  —  1)  to  indicate  the  set  1, 2, . . . .  <  -  1 . 

The  assumption  that  lAil,|A2l  <  1  results  in  the  distribution  of  ^  ^  ^  being  stationary 
and  it  is  the  distribution  unconditional  on  other  plots  that  we  use  when  at  the  beginning 
of  a  block.  We  are  assuming  here  that  while  the  plot  at  the  beginning  of  a  new  block  is 
not  adjacent  in  the  field  to  the  previous  plot,  the  two  plots  are  relatively  close.  Thus  while 
the  two  plots  have  similar  characteristics  (i.e.  the  sime  marginal  distribution)  they  can  be 
assumed  independent. 

The  intuition  of  a  smoothly  varying  fertility  suggests  that  both  A,  and  A2  should  be 
non-negative.  The  restriction  Ai  >  A2  is  enforced  to  aid  identifiability  as  ’h.cr  marginal 
distribution  of  Yt  remains  unchanged  if  we  interchange  Aj  and  A2.  Inde'^j,  the  marginal 
distribution  of  the  fertility  level  Ft  also  remains  unchanged,  with  this  interchange  affecting 
only  the  distribution  of  the  fertility  gradient  GV  The  choice  of  Aj  being  tlie  laiger  coincides 
with  our  intuition  governing  the  definition  of  F  and  G. 

In  line  with  e  representing  measurement  error,  we  shall  assume  that  the  Ct,  {t  =  1 , . . . ,  noplo 
are  independent  (and  independent  of  F),  and  further  make  the  distributional  assumption. 


e,  ~  t  =  I,. noplots 


>  0. 
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While  e  could  be  absorbed  into  F,  we  make  the  distinction  as  the  term  e  is  directly 

interpretable  as  measurement  error  and  observed  data  often  supports  its  presence  (Whittle. 

1954).  The  distribution  of  y  is  therefore  dependent  on  the  parameter  Aj,  A2), 

or  e  =  (o■^%/er^I,  A,,  A2)  where  a'^  =  A;)(i-A,^)‘^grad  +  's  the  variance  of  any 
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observation  y  given  D  and  r  only,  and  %fert  =  100  (1  -^)  >s  the  percentage  of  this 
variance  explained  by  the  fertility. 

The  condition  0  <  Aj  <  Aj  <  1  forms  a  triangle  in  (Ai,A2)  parameter  space,  with  an 
open  side  (being  the  bound  to  non-stationarity)  and  two  closed  sides.  This  suggests  two 
natural  special  cases  for  modelling  field  fertility,  with  the  more  general  parameter  values 
lying  between  these  special  cases.  (We  ignore  the  non-stationary  side;  see  Section  6).  The 
first  special  case  is  obtained  by  forcing  A2  =  0.  This  model,  while  assuming  that  the  fertility 
levels  at  consecutive  plots  are  related,  assumes  that  the  fertility  gradients  are  independent. 
This  typically  results  in  jagged  fertility  trends.  The  second  case  results  by  forcing  A,  =  A2 
and  potentially  produces  smoother  fertility  trends  than  the  first. 

Figure  1  displays  these  two  special  cases  in  relation  to  our  (Aj,  A2)  parameter  space,  to¬ 
gether  with  various  fertility  models  previously  proposed  for  field  fertility.  These  include  the 
simplest  conditional  auto- regression  (Besag,  1974),  the  simplest  simultaneous  auto- regression 
(Whittle,  1954),  the  first  difference  model  (Besag  and  Kempton.  1980),  the  second  differ¬ 
ence  model  (Green,  Jennison  and  Seheult,  1985),  and  the  ARIM.A  (1.1.0)  (Glceson  and 
Cullis,  1987),  AR(1)  (Patterson,  1983)  and  white  noise  models. 

Figure  1;  ***  Figure  1  about  here 

Our  fertility  model  is  a  special  case  of  the  AR(2)  model,  with  A,  and  A2  being  the  roots 
of  the  characteristic  equation.  Forcing  these  roots  to  be  real  and  positive  eliminates  models 
with  oscillating  correlograms  that  can  take  negative  values.  We  argue  that  these  phenomena 
do  not  conform  with  the  notion  of  a  fertility  trend,  where  the  correlations  sltould  decrease 
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with  increasing  distance  and  remain  non-negati'’e.  It  present,  it  is  assumed  that  they  can  be 
more  directly  modelled  (such  as  with  direction  of  harvester  as  a  covariate),  and  not  included 
with  the  fertility. 

The  model  with  Aj  =  is  the  AR(2)  model  with  equal  roots  of  the  characteristic 
equation.  The  resulting  autocorrelation  function  has  the  form  pk  =  (l  +  ^i-  This 

autocorrelation  function  is  unlike  those  for  any  other  AR(2)  model  in  that  while  it  is  monot- 
ically  decreasing  away  from  zero,  it  is  flat  at  the  origin  rather  than  spiked.  This  special 
case  of  the  AR(2)  model  was  not  studied  by  Box  and  Jenkins  (1976),  for  example.  The 
same  autocorrelation  structure  also  arises  from  what  may  be  considered  the  simplest  truly 
bilateral  simultaneous  autoregression,  namely  Ft  =  a(F,_i  +  Ft+\)  +  where  a  — 
Whittle  (1954)  considered  the  two-dimensional  version  of  this  process  and  described  data 
that  supported  this  phenomenon. 

2.2  Bayesian  Estimation 

To  calculate  the  posterior  distribution  of  r  we  must  calculate  the  integral 

p{r\Y)  =  jv{T,e\Y)dO 

where  p(K|r,^)  is  the  likelihood  and  p{t,9)  =  p{t\0)p{6)  is  our  prior  for  (r,  (1).  In  general 
this  integral  cannot  be  evaluated  analytically,  and  so  a  numerical  technique  must  be  used. 
The  computation  can  be  simplified  by  using  a  multivariate  normal  prior  for  {r\0)  (or  with 
simple  modifications  a  mixture  of  multivariate  normals).  It  then  follow.s  that  t?)  oc 

p(  V’lr,  0)p{t\6)  is  also  a  multivariate  normal,  and  p(T|y')  =  / p{t\Y,6)p(0\Y )  dO  is  a  mixture 
of  the.sc  normal  distributions. 

The  mean  and  variance  of  {t\Y,6]  and  the  likelihood  given  0.  p(}  Id),  can  l)e  calculated 
directly  using  the  Kalman  filter  with  slate  A'(  =  (F(,G'(."i . 'nmiMis)  ♦’ven  in  the  presence 
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of  missing  values  (Kalman,  1960).  This  suggests  the  use  of  importance  sampling  to  estimate 
the  posterior  for  t  by  first  simulating  {i=l, . . . ,  N)  independently  from  the  sampling 
importance  density  f{0),  and  then  forming  the  estimate 


(4) 


The  distribution  f{6)  is  chosen  such  that  it  can  be  easily  simulated  from,  and  is  as  close  to 
the  unknown  p{9\Y)  as  possible.  As  N  approaches  infinity,  p(T|y')  approaches  p{t\Y). 
Another  approximation  is 

PiT\y)  =  PiT\6,Y),  (5) 


where  6  is  the  value  of  6  that  maximises  the  posterior  distribution  of  p(^|V’).  A  numerical 
maximisation  of  p{6\Y)  is  therefore  required  to  obtain  6. 

Figure  2  displays  the  posterior  of  a  treatment  contrast  (early  spraying  —  no  spraying) 
in  the  mildew  control  trial  (Draper  and  Guttman,  1985;  Jenkyn  et.  ai,  1979)  estimated  by 
equations  (4)  and  (5)  under  the  prior  p{t,6)  oc  Similar  results  were  obtained  with  other 
vague  priors.  While  conditioning  on  6  resulted  in  a  good  approximation  with  considerably 
less  computation,  as  expected  it  underestimates  the  uncertainty  about  the  value  of  this 
contrast.  It  obtains  the  same  mean  but  its  variance  is  underestimated  by  about  19%. 


Figure  2:  ***  Figure  2  about  here  *** 


The  accuracy  of  the  approximation  (5)  is  not  surprising  when  one  notes  its  similarities 
with  Restricted  Maximum  Likelihood  (REML)  estimation  (Patterson  and  Thompson,  1971; 
Laird  and  Ware,  1982).  When  the  prior  p(t, 0)  oc  1  is  used,  0  equals  the  REML  estimate  of 
0  and  the  (1— a)  highest  posterior  density  region  of  the  approximate  posterior  distribution 
of  T  equals  the  (1— a)  confidence  interval  from  conditioning  on  the  REML  estimate  of  0. 
To  correct  the  underestimation  of  uncertainty  that  occurs  in  REML,  it  has  been  suggested 
that  the  REML  estimate  of  be  inflated  by  ^  where  n  is  the  number  of  observations  and 
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df =n  —  d\m{0) .  However,  here  this  increcises  the  variance  by  only  about  9%  compared  with 
the  19%  that  would  be  required  to  make  it  exact. 

Note  that  a  better  choice  of  prior  than  p{t,6)  a  1  is  likely  to  be  available.  Even  ig¬ 
noring  the  fact  that  6  contains  variance  terms  for  which  a  non-uniform  prior  may  be  more 
appropriate,  for  many  trials  prior  information  about  r  will  be  available.  This  is  particularly 
likely  to  be  the  case  when  r  represents  the  effect  of  different  varieties.  In  early  generation 
variety  trials  the  prior  for  r  should  reflect  the  genetic  relationship  between  the  varieties  (see 
for  example  Cullis  et.  al.,  1990).  For  later  trials,  the  information  derived  from  earlier  trials 
with  these  varieties  can  be  utilised. 

An  alternative  to  using  a  pre-specified  prior  is  to  use  a  random  effects  model  where  the 
variety  effects  have  an  exchangeable  joint  distribution.  The  prior  distribution  of  the  variety 
effects  can  be  estimated  from  the  data,  leading  to  an  empirical  Bayes  approach.  If  that 
distribution  has  a  parametric  form,  we  have  a  parametric  empirical  Bayes  model  (Morris, 
1983).  A  simple  parameterisation  of  the  distribution  for  this  population  is  to  assume  t 
multivariate  normal  with  mean  /Xt-I  and  variance  matrix  a^Cr-  Here  Cr  is  a  pre-specified 
covariance  matrix  for  r,  1  is  the  notreats  vector  of  ones,  and  /Xt  and  are  unknown  scalars 
that  could,  for  example,  be  estimated  by  the  mode  of  their  posterior  density. 

Figure  3  illustrates  the  effect  of  using  a  random  effects  model  instead  of  a  fixed  effect 
model  on  the  ARC  2  trial  (a  trial  with  38  treatments  replicated  3  times  and  nothing  unusual 
in  the  way  of  outliers  or  fertility  jumps).  Here  we  have  used  the  approximation  just  described 
with  Ct  equal  to  the  identity  and  the  model  with  Xi  =  X2-  Note  that  the  random  effects  model 
results  in  estimated  treatment  effects  (posterior  means)  with  a  smaller  range.  The  standard 
deviations  of  the  treatment  difference  posteriors  were  similarly  reduced.  This  shrinkage 
agrees  with  our  intuition  that  the  treatment  with  the  largest  estimated  effect  from  the  fixed 
effect  model  is  likely  to  be  overestimated.  The  reason  for  this  treatment  's  favourable  result 
is  at  least  in  part  likely  to  be  due  to  random  deviations  (as  it  would  be  completely  if  the 
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true  variety  effects  were  equal). 


Figure  3:  ***  figure  3  about  here*** 


3  Modelling  Outliers  and  Fertility  Jumps 

3.1  The  Robust  Model 

We  now  extend  our  model  to  allow  for  the  possibility  of  outliers  and  fertility  jumps.  The 
term  outlier  is  used  to  refer  to  a  single  observation  whose  value  was  generated  by  a  different 
mechanism,  while  a  fertility  jump  refers  to  a  sudden,  relatively  large  shift  in  the  fertility 
level  compared  to  the  majority  of  the  observations.  Both  occurrences  are  considered  rare, 
but  when  present  can  greatly  influence  the  resulting  inferences.  We  retain  the  decomposition 
in  equation  (1),  and  the  treatments  r  and  design  matrix  D,  but  assume  that  a  fertility  jump 
modifies  the  distribution  of  F  and  an  outlier  modifies  the  distribution  of  e. 

The  field  fertility  is  modelled  as  in  equations  (2)  and  (3).  However,  given  that  there 
is  a  fertility  jump  at  plot  t,  we  assume  that  the  fertility  at  plot  t  is  independent  of  the 
fertility  at  plot  t  —  1.  We  apply  the  same  justification  as  was  used  in  Section  2  when  plot  t 
is  at  the  beginning  of  a  block,  and  hence  apply  the  same  distributional  assumptions.  Thus 
equation  (3)  is  used  if  plot  t  is  either  at  the  beginning  of  a  block  or  at  a  fertility  junip.  with 
equation  (2)  being  used  otherwise. 

We  model  outliers  by  modifying  the  distribution  of  the  measurement  errors  t.  While 
retaining  normality  and  independence  of  the  Cj,  we  inflate  the  variance  of  t(  by  a  constant 
factor,  when  observation  i  is  an  outlier.  We  shall  henceforth  assume  this  variance 
inflation  factor  to  be  100  (see  Section  5  for  a  discussion  of  this  value).  Thus,  while  an  outliei 
observation  y,  still  has  a  distribution  centred  on  (Z^t),  +  its  distribution  is  considerably 
more  spread  out. 
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To  avoid  a  possible  identiflability  problem,  we  assume  that  a  fertility  jump  cannot  occur 
at  the  beginning  of  a  block.  We  also  make  the  assumption  that  a  fertility  jump  and  an 
outlier  cannot  both  be  present  at  the  same  plot.  This  simplifies  our  model  by  allowing  at 
most  three  possible  cases  for  each  plot:  a  fertility  jump,  an  outlier,  or  neither.  The  fourth 
possibility  of  both  a  fertility  jump  and  an  outlier  should  be  uncommon,  and  we  note  in  this 
case  that  the  distribution  of  y^  is  affected  little  by  the  exact  value  of  the  fertility  at  plot  t. 
Hence  in  this  case  we  re-define  the  fertility  jump  at  plot  t  to  be  at  plot  t  +  \  .  While  this 
results  in  an  incorrect  estimate  of  the  plot  fertility  at  plot  f,  we  note  this  and  consider  it  to 
be  of  little  consequence,  as  the  posterior  distribution  of  r  which  is  of  primary  interest  will 
be  affected  negligibly. 

Formally,  we  shall  let  Ct  denote  the  condition  of  plot  t.  Thus  Cj  =  1  if  plot  i  has  neither 
an  outlier  nor  a  fertility  jump,  C(  =  2  if  plot  t  has  an  outlier,  and  c,  =  3  if  plot  f  is  at  a 
fertility  jump.  Then  c  =  (cj, . . . ,  c„op/otj)  denotes  the  condition  of  all  the  plots.  We  define 
the  submodel  to  be  the  model  with  plot  conditions  c. 

Our  model  for  Y  consists  of  a  mixture  of  these  submodels  where  each  submodel  is 
weighted  according  to  the  prior  probability  of  model  Me-  The  likelihood  for  Y  is  therefore 

piY)  = 

C 

and  the  posterior  distribution  of  r,  the  quantity  of  interest,  is  given  by 

p(r\Y)cxY,P{r\Y.M,)p(MAY). 

C 

We  assume  that  the  prior  specifies  the  conditions  of  the  different  plots  to  be  statistically 
independent  and  that  the  conditions  1,2,  and  3  have  prior  probabilities  96%,  2%,  and  2% 
respectively.  These  choices  are  discussed  in  Section  5. 
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3.2  Comparison  with  robust  time  series 

In  the  absence  of  the  term  -Dt,  we  obtain  a  robust  model  similar  to  those  employed  in 
robust  non-seasonal  time  series  (Harrison  and  Stevens,  1976;  West  and  Harrison,  1990). 
Estimation  under  such  models  is  aided  by  the  fact  that  an  observation  is  most  correlated 
with  observations  on  nearby  plots,  the  neighbours.  Thus  the  condition  of  a  ])lot  can  be  well 
estimated  with  knowledge  only  of  the  observations  on  the  neighbouring  plots.  Estimation 
can  then  proceed  recursively  by  updating  the  plots  in  field  order,  so  that  a  plot  is  considered 
at  about  the  same  time  as  its  neighbours. 

In  the  presence  of  the  term  Dr,  the  neighbourhood  structure  described  above  is  more 
complex,  with  the  neighbouring  plots  not  only  being  those  nearby  in  the  field,  but  also  those 
receiving  a  similar  treatment  combination.  A  recursive  procedure  is  therefore  not  as  efficient, 
as  it  becomes  difficult  to  update  the  plots  in  such  a  way  that  a  plot  is  updated  at  about  the 
same  time  as  its  neighbours.  The  simplest  solution  to  this  problem  would  be  to  condition  on 
an  estimate  of  r,  and  then  estimate  the  condition  of  the  plots  as  before.  Figure  4  compares 
the  results  of  using  maximum  likelihood  on  the  full  model  to  estimate  r  with  the  exact 
Bayesian  solution  in  some  simple  cases. 

Figure  4:  ***  figure  4  about  here  *** 

In  Figure  4  we  assumed  that  F=0  and  ,  so  that  the  observations  are  independent 

with  a  variance  of  1,  or  of  100  if  conditioned  to  be  an  outlier.  The  (scalar)  mean  r  was 
assumed  to  have  the  improper  uniform  prior  and  only  ,.he  first  observation  had  positive  prior 
probability  (2%)  of  being  an  outlier.  The  approximation  is  good  with  nine  replicates,  but 
seriously  underestimates  the  uncertainty  about  whether  the  observation  is  an  outlier  or  not 
with  three  replicates.  The  problem  occurs  because  the  likelihood  is  the  mixture  of  the  two 
normal  likelihoods  corresponding  to  the  two  submodels,  and  with  few  replicates  the  critical 
factoi'  is  which  of  these  likelihoods  dominates  the  other.  Maximum  likelihood  cannot  be 
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expected  to  perform  well  when  the  log-likelihood  is  so  noii-qu^idratic.  nor  is  it  tri\ial  to 
ensure  that  the  global  maximum  and  not  just  a  local  maximum  is  achieved. 

3.3  Estimation  for  the  Robust  Model 

The  over-simplified  model  used  for  Figure  4  suggests  that  the  posterior  probability  of  an 
observation  being  an  outlier  is  well  approximated  by  zero  unless  its  cross-validated  residual 
is  at  least  3.  Since  such  a  large  residual  is  rare  under  the  non-robust  model,  this  suggests  that 
most  of  the  submodels  of  our  model  will  have  very  low  posterior  probability.  We  therefore 
propose  to  calculate  only  the  submodels  that  have  non-negligible  posterior  probabilities 
M<.,c  G  C,  and  approximate  the  posterior  probabilities  for  the  other  submodels  by  zero.  We 
defer  the  estimation  of  C  to  Section  3.4,  but  now  take  it  as  given. 

We  approximate  the  posterior  distribution  of  r  as  follows; 

V{r\Y)  «  Y.lv{r.B,M,\Y)dO 

cec-' 

«  Epi^c\Y)p(T\LMc,Y).  (G) 

cCC 

In  equation  (6),  0^  is  the  value  of  6  that  maximises  p{9\Mc,Y),  and 

p{M,\Y)  cx  p(r|M,)p(A/,) 

-  p{M,)Jp{Y\0,MM0\M,)d0 

^  p{mmy\Lmm0c\m,)^  {-) 

where  p{0\Mc)  is  the  prior  for  0  under  submodel  Me,  and  is  typically  assumed  to  be  the 
same  for  every  submodel.  Thus  each  submodel  is  estimated  separately,  as  in  Section  2.J. 
and  the  icsults  combined  under  the  assumption  that  integrating  0  out  of  p{)'.0\M^)  can  be 
well  airproximated  by  maximising  it  with  respect  to  0. 

The  simpler  approximation  that  conditions  on  a  single  value  0  of  0  that  maximisc's 
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E,6CP(«,W,|i'), 

Pir\Y)  =  fZPiY«,Kir)M 

C 

«  / Y.p{A&^M,,Y)p{M^d^-)p[e\Y)de 

c6C 

Y,p{t\6,M,,Y)p[MA6,Y)  (8) 

cec 

applies  the  same  approximation  as  in  Section  2.2  to  the  complete  model.  A  comparison  of 
these  two  approximations  is  illustrated  in  Figure  5,  where  we  again  assume  the  simplified 
model  with  zero  fertility,  so  that  and  the  observations  are  independent.  Observations 

receiving  the  first  treatment  were  all  assumed  to  be  zero  except  for  i/],  while  observations  on 
the  other  treatments  were  such  that  their  residual  sum  of  squares  equalled  their  degrees  of 
freedom  (so 

Figure  5:  ***  Figure  5  about  here  *** 

In  this  case  approximation  (7)  is  very  good  (and  is  exact  if  the  prior  p{(^oh%)  oc  1  is 
used  in  (7)  in  place  of  the  assumed  p{(T^\^)  oc  l/cTp^jj).  Approximation  (8)  underestimates 
the  uncertaint}'  about  whether  the  first  plot  contains  an  outlier  or  not,  especially  for  small 
trials. 

3.4  Determining  the  Conditioning  Set  of  Outliers  and  Jumps 

Due  to  the  large  number  of  possible  submodels,  it  is  likely  that  any  practical  choice  of  C. 
the  conditioning  set  of  anomalies,  will  in  total  comprise  a  small  percentage  of  the  posterior 
probability;  the  majority  of  the  probability  being  spread  very  thinly  over  the  \er\-  large 
number  of  remaining  submodels.  Nevertheless,  the  posterior  for  the  treatment  difference 
T.  —  Tj  will  primarily  be  influenced  by  the  presence  of  outliers  on  any  plots  with  treatments 
i  or  j  applied,  or  by  fertility  jumps  on  plots  close  to  these  plots.  It  is  the  marginal  posterior 
probabilities  of  outliers  and  fertility  jumps  on  these  plots  that  is  of  primary  importance 
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when  calculating  the  posterior  for  this  treatment  difference.  These  marginal  posteriors  can 
be  reasonably  approximated  by  a  few  submodels  since  to  a  first  approximation  the  posterior 
occurrences  of  outliers  and  fertility  jumps  can  be  assumed  independent.  This  approximation 
is  most  likely  to  be  violated  when  the  addition  of  an  outlier  or  fertility  jump  greatly  modifies 
the  estimate  of  influencing  the  probabilities  of  future  outliers  or  fertility  jumps  through 
different  distributional  assumptions  (likely  only  in  small  trials  or  when  the  original  outlier 
or  fertility  jump  has  very  high  posterior  prebability),  or  when  the  outliers  or  fertility  jumps 
in  question  are  close  under  the  neighbourhood  structure. 

To  determine  the  set  of  submodels  Me,  c  E  C,  that  have  non-negligible  posterior  proba¬ 
bility  we  adopt  a  recursive  procedure  beginning  with  the  submodel  Ct=l  Vt.  This  submodel 
corresponds  to  no  outliers  or  fertility  jumps  and  often  has  high  posterior  probability;  we 
adopt  the  convention  of  always  including  it  in  C. 

Given  that  we  have  just  included  the  submodel  Me  in  C,  we  then  consider  the  submodels 
Me'  in  Ce'  which  have  the  same  outliers  and  fertility  jumps  as  Me  plus  one  extra  outlier 
or  fertility  jump.  The  submodels  in  Ce'  are  then  approximately  ranked  according  to  their 
posterior  probability,  and  estimated  in  turn  until  a  submodel  from  Ce'  is  rejected  from  C  due 
to  its  low  posterior  probability.  Here  the  submodel  Me’  is  rejected  if  its  posterior  probability 
is  less  than  a  proportion  q  of  the  sum  of  the  posterior  probabilities  for  the  submodels 
previously  considered  from  C'  and  Me-  The  proportion  q  is  taken  to  be  small,  and  in  our 
examples  we  took  a  =  2%.  Choosing  a  to  be  close  to  the  (usually  small)  prior  probability  of 
an  outlier  or  a  fertility  jump  seems  reasonable  because  it  implies  roughly  that  when  the  data 
alone  provide  evidence  against  a  plot  being  an  outlier  or  a  jump,  that  possibility  is  ignored. 
(Of  course,  the  posterior  probability  of  an  outlier  or  a  jump  is  reduced  by  also  taking  account 
of  the  prior  knowledge  that  these  are,  by  definition,  relatively  rare  events).  The  recursion  is 
continued  by  considering  those  submodels  with  an  extra  outlier  or  fertility  jump  from  those 
submodels  M<->  accepted  into  C. 
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There  are  various  ways  of  determining  the  ranking  of  Likely  outliers  can  be  detected 
by  the  estimated  measurement  errors  with  large  magnitude  from  the  submodel  (note 
that  Oc'  is  likely  to  be  close  to  Oc  due  to  the  similarity  of  the  submodels).  Similarly.  likel\ 
fertility  jumps  may  be  detected  by  values  of  — +  from  the  submodel  A 
more  refined  estimate  for  a  selection  of  submodels  could  be  baised  on  V)  where 

fc  is  the  estimate  of  r  from  M^,  or  even  better  based  on  p{Mc'\6c,y)-  These  rankings  for 
Cc'  are  likely  to  be  progressively  more  accurate  but  require  progressively  more  computation. 
They  all  require  considerably  less  computation  than  our  final  estimate,  p(Mc')0cS  1  ),  which 
requires  a  numerical  optimisation  and  should  only  be  required  for  a  small  subset  of  submodels 
in  Cc'- 

4  Examples 

We  now  consider  three  trials  as  examples.  All  are  early  grain  variety  trials  conducted  by 
the  Agricultural  Research  Council  of  Great  Britain  and  consist  of  3  or  4  blocks  with,  each 
variety  (treatment)  applied  to  exactly  one  plot  in  each  block.  The  first  two  trials  (.ARC  S 
and  ARC  6)  illustrate  respectively  the  case  where  there  is  uncertainty  about  whether  a  single 
outlier  or  fertility  jump  is  present,  while  the  third  trial  (ARC  1)  illustrates  the  occurrence 
of  many  possible  outliers  and  fertility  jumps. 

4.1  A  Possible  Outlier:  Trial  ARC  8 

The  model  with  general  (Ai,A2)  provided  negligible  gain  in  likelihood  p{Y\d)  over  either  of 
the  two  special  crises  A2=0  or  Ai  =  A2,  and  we  therefore  restrict  consideration  to  the  latter  two 
models.  Figure  6  displays  the  treatment  adjusted  data  (i.e.  Y  —  Dt  where  f  is  the  estimate  of 
T  under  the  submodel  C(  =  l  Vi)  in  field  order  for  the  three  blocks  with  the  treatment  number 
as  the  plotting  symbol.  The  lines  are  the  estimated  fertilities  from  the  A,  =  A2  model  (solid) 
and  the  A2=0  model  (dotted).  While  the  posteriors  for  r  were  very  similar  under  the  two 


15 


-V 


models,  the  estimated  fertility  is  smoother  under  the  ^,=^2  model,  with  correspondingly 
larger  estimated  measurement  errors. 

Figure  6:  ***  Figure  6  about  here  *** 

Figure  7:  ***  Figure  7  about  here  *** 

The  large  magnitude  of  the  estimated  measurement  error  on  plot  23  suggests  that  this 
value  is  atypical,  and  Figure  7  is  similar  to  Figure  6  but  assumes  the  submodel  C23=2.  c,  =  l 
Vt^23.  This  is  almost  equivalent  to  removing  this  observation  from  the  analysis.  Table  1 
summarises  the  estimation  of  these  two  submodels  under  the  two  models. 

Table  1:  ***  table  1  about  here  *** 

Both  models  achieve  similar  loglikelihoods,  and  this  together  with  the  fact  that  for  either 
submodel  the  models  achieve  similar  posteriors  for  r  suggests  that  for  this  trial  submodel  un¬ 
certainty  is  more  important  than  model  uncertainty.  With  a  prior  probability  of  an  outlier  of 
2%,  this  translates  (by  approximation  (7))  into  posterior  probabilities  of  an  outlier  on  plot  23 
of  82%  and  79%  for  the  ^2=0  and  Ai=A2  models  respectively.  Approximation  (S)  how’ever 
results  respectively  in  posterior  probabilities  of  98%  and  97%,  and  essentially  conditions  on 
plot  23  containing  an  outlier. 

For  comparison,  a  complete  Bayesian  analysis  with  prior  p{d)  oc  1/cr^  resulted  in  a  pos¬ 
terior  probability  of  C23=2  of  about  68±4%  (calculated  using  importance  sampling  as  in 
Section  2.2).  Thus  approximation  (7)  also  underestimates  the  uncertainty  about  the  condi¬ 
tion  of  plot  23,  but  not  as  much  as  approximation  (8). 

Finally  we  note  that  it  is  the  possibility  of  C23=l  or  2  that  is  the  primary  issue  for  this 
trial.  The  next  most  likely  outlier  appears  to  be  on  plot  33  with  a  posterior  probability  of 
less  than  3%,  and  the  most  likely  fertility  jump  probably  occurs  at  plot  43,  wdth  a  posterior 
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probability  of  about  6%.  These  possibilities  have  little  influence  on  the  posterior  for  t.  The 
effect  on  the  posterior  of  t  of  uncertainty  about  the  condition  of  plot  23  will  be  illustrated 
in  Section  5. 

4.2  A  Possible  Fertility  Jump:  Trial  ARC  6 

Whereas  the  trial  ARC  8  was  dominated  by  a  large  measurement  error  relative  to  the  fertility 
trend,  the  trial  ARC  6  provides  the  contrasting  situation  where  is  very  small.  Here  we 
shall  assume  that  =  0  for  trial  ARC  6  to  illustrate  a  possible  fertility  jump.  Figure  8 
graphs  the  treatment  adjusted  yield  (under  the  model  with  A]=A2)  against  plot  number  in 
the  same  manner  as  for  trial  ARC  8.  Note  that  here  there  are  four  blocks,  and  since  =  0 
the  estimated  fertility  equals  the  treatment  adjusted  yield. 

Figure  8;  ***  Figure  8  about  here  *** 

Table  2  summarises  the  estimation  of  the  two  submodels  C(  =  l  V<  and  C5=3,  C(  =  l 
under  the  two  models  with  cr^bs^O  (note  that  only  when  Ai  =  A2  and  C5=3  would  a^bs 
estimated  to  be  non-zero). 

Table  2:  ***  table  2  about  here  *** 

Unlike  in  the  previous  example,  the  choice  of  model  now  influences  the  evidence  for  the 
submodels;  the  fertility  jump  at  plot  5  is  more  likely  under  the  model  with  A,  =  A2.  Assuming 
a  prior  probability  of  a  fertility  jump  of  2%,  the  posterior  probability  that  C5=3  is  63%  when 
Ai=A2  compared  to  16%  when  A2=0.  This  increcised  model  influence  is  to  be  expected  due 
to  the  lack  of  measurement  error,  making  the  distribution  of  V  more  dependent  on  the 
distribution  of  F . 

We  note  that  while  the  model  A2=0  was  0.6  units  of  loglikelihood  superior  to  the  model 
Ai=A2  when  conditioning  on  C5  =  l,  it  is  0.2  units  inferior  when  taken  unconditionally.  Thus, 
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when  the  possibility  of  fertility  jumps  is  introduced,  there  is  quite  a  shift  in  evidence  about 
which  fertility  model  is  appropriate.  Furthermore,  under  the  submodel  C(  =  l  V/,  it  is  the 
posterior  for  the  treatment  difference  T17  —  rjg  that  differs  most  between  the  two  models.  It 
is  also  this  treatment  difference  that  is  most  influenced  b\'  the  introduction  of  the  fertility 
jump  at  plot  5  (since  treatments  17  and  16  were  applied  to  plots  5  and  4  respectively), 
and  we  shall  show  in  Section  5  that  allowing  for  the  possibility  of  fertility  jumps  makes  the 
treatment  posteriors  under  the  different  models  more  similar. 

4.3  Multiple  Possible  Outliers  and  Fertility  Jumps:  Trial  ARC  1 

The  treatment  adjusted  data  is  plotted  in  Figure  9  together  with  the  estimated  fertilit}- 
from  the  Ai  =  A2  model.  A  dominant  feature  here  is  the  inability  of  the  estimated  fertility 
to  increase  as  quickly  as  the  treatment-adjusted  data  over  plots  1  to  10.  A  fertility  jump 
at  plot  5  appears  more  appropriate,  and  Figure  10  displays  the  results  conditional  on  the 
submodel  Cs=3,  C(  =  l  This  fertility  jump  is  well  supported  by  the  data,  with  a  posterior 

probability  of  over  95%  from  a  prior  of  only  2%. 

Figure  9;  ***  Figure  9  about  here  *** 

Figure  10:  ***  Figure  10  about  here  *** 

We  shall  restrict  attention  here  to  the  smoother  model  Xi=X2  because  while  it  is  superior 
to  the  A2=0  model  by  only  0.1  units  of  loglikelihood  when  conditioned  on  Cs^l,  it  is  superior 
by  1.4  units  when  C5=3.  Due  to  the  dominance  of  05=3  over  C5  =  l  (under  either  model),  the 
data  strongly  favours  the  Ai  =  A2  model  under  the  robust  model. 

Table  3  shows  the  13  submodels  with  highest  approximate  posterior  probability,  and  their 
posterior  probabilities  conditional  on  C  only  containing  these  submodels.  Table  4  shows  the 
marginal  posterior  probabilities  of  the  7  most  likely  outliers  and  fertility  jiim|xs  from  our 
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analysis  (with  q=2%,  see  Section  3.4).  Note  that  conditioning  on  a  single  submodel  is  not 
appropriate,  a5  the  submodel  with  highest  posterior  probability  conditions  on  a  fertility  jump 


at  plot  5  only,  while  the  marginal  posterior  probability  of  an  outlier  on  plot  49  is  63%. 

Table  3:  table  3  about  here  *** 

Table  4;  *'*  table  4  about  here  *** 


5  Sensitivity  to  Model  Specification 

Due  to  the  large  number  of  submodels,  it  is  not  practical  to  compute  the  likelihood  as  a 
function  of  the  parameters  pom  and  Pjump-  Furthermore,  the  likelihood  function  can  be  quite 
flat  indicating  that  the  data  contain  little  information  concerning  these  parameters,  and 
we  have  conditioned  on  values  of  these  parameters  (2%)  chosen  prior  to  data  collection. 
Nevertheless,  the  value  of  these  parameters  can  have  a  substantial  effect  on  the  posterior  for 
T,  especially  if  they  are  set  too  low  in  the  presence  of  outliers  and  fertility  jumps. 

It  is  therefore  useful  to  examine  how  the  posterior  for  r  is  influenced  by  our  prior  choice 
for  Pout  and  pjump,  and  this  information  is  easily  available  under  the  approximation  (6) 
and  (7).  We  note  that  under  this  approximation  the  submodel  priors  influence  the  posterior 
for  T  only  through  the  mixing  proportions  of  the  posteriors  for  t  based  on  each  submodel 
(assuming  the  same  C).  Figure  11  displays  the  posterior  for  the  treatment  differences  Ts  -  t- 
and  Tie  ~  Ds  in  trials  8  and  6  as  a  function  of  pout  and  Pj„mp  respectively.  The  left  graph  in 
Figure  11  assumes  that  C  contains  only  the  two  submodels  c,  =  l  V<  and  C23=2,  C(  =  l,  yt^'23 
and  the  model  A2=0.  The  other  two  graphs  in  Figures  11  both  assume  that  C  contains  only 
the  submodels  C(  =  l  Vi  and  C5=3,  C(  =  l  Vi/5.  The  center  graph  in  Figure  11  assumes  the 
model  Ai  =  A2  while  the  model  A2=0  is  assumed  on  the  right.  In  all  cases  the  ^%,  2^%,  97^%. 
and  99^%  quantiles  are  graphed  along  with  the  posterior  mean. 
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Figure  11:  Figure  11  about  here 
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The  major  feature  of  these  graphs  is  the  insensitivity  of  the  posterior  to  the  precise  value 
of  Pout  or  Pjump-  For  these  examples,  any  value  of  Pout  or  Pjump  between  and  5%  will  yield 
relatively  similar  posteriors  compared  to  the  extremes  of  0%  and  100%  (which  correspond  to 
the  two  submodels).  The  two  graphs  to  the  right  in  Figure  11  also  illustrate  how  similar  the 
posteriors  from  different  models  can  be  after  taking  into  account  the  submodel  uncertainty. 
The  indicated  quantiles  are  similar  when  Pjump=2%  despite  their  being  quite  different  for 
each  submodel. 

We  now  examine  the  influence  of  (used  to  model  outliers)  on  the  posterior  proba¬ 
bilities  of  the  submodels.  As  above  for  poui,  while  it  is  intended  that  a  value  of  kom  fixed 
in  advance  be  used,  it  is  useful  to  see  how  the  posterior  for  r  varies  aus  a  function  of  ko,„. 
Unfortunately  the  situation  is  not  so  simple  when  varying  k^uiy  as  each  submodel  must  be 
re-fitted  (i.e.  optimised  over  0)  with  every  new  value  of  /.'out-  In  general  this  will  require 
considerable  computation,  so  instead  we  provide  a  simple  approximation. 

By  definition  of  an  outlier  we  require  /:out>l,  and  in  practice  we  expect  it  to  be  at  least  as 
large  as  5  (which  will  approximately  remove  an  outlier  from  the  analysis).  If  we  approximate 
the  likelihood  by  conditioning  on  estimates  of  t  and  0,  then  the  ratio  of  the  likelihoods  for 
Me  under  /:out  =  /^outi  to  koua  is  approximately. 
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-'"oul2‘^obs.2  1=1 


1=1  / 


(9) 


\  fcoull  d'obs.l  y  V  ^/.'outl^^obs.l  1=1 

where  Me  conditions  plots  o(l),o(2), . . .  ,o(n)  only  as  outliers,  is  the  estimated  mea¬ 

surement  error  on  plot  o(z)  under  /:oui  =  ^’ouij7  and  is  the  estimate  of  when  /'oui  =  /'..uij 
This  approximation  is  reasonable  if  /'ouiC^obs  ‘s  large  compared  with  (so  that  outliers  are 
approximately  independent  of  the  other  plots),  which  will  typically  be  the  case. 

Expression  (9)  still  depends  on  the  and  to  avoid  the  estimation  of  tliesc  for  every 

value  of  fcouti  we  could  condition  on  the  same  value  of  CTobs  obtained  from,  say,  /:o„i  =  10.  The 
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resulting  approximation, 


2JL2  -2  IZ^i.o(0 ) 

^^outl^obs  1=1  / 


A  -2 

9p  -  2  2^  ‘■2,o(.) 

“^OUt2^ob5  1=1 


is  very  good  and  easily  used  to  evaluate  the  effect  of  changing  k^ui-  Note  that  this  assumption 
of  similar  6  is  more  reasonable  here  than  when  comparing  different  submodels  as  practical 
values  of  fcom  will  all  essentially  remove  outliers  from  the  analysis.  Furthermore,  for  values 
of  |£_,,o(i)|  of  about  4  or  5  that  create  submodel  uncertainty,  the  ratio  of  exponentials  in 
equation  (10)  is  approximately  unity  (especially  for  larger  A'oui),  which  results  in  the  ratio  of 
likelihoods  being  approximately  ^out2/^outi- 

Thus  this  approximation  suggests  using  the  (pouu^out)  combination  of  (pout2|‘^i  ^'oun ) 
instead  of  (pout2,^out2)i  and  so  Figure  11  can  again  be  used.  This  approximation  also  illus¬ 
trates  how,  for  larger  /cout,  it  is  the  ratio  of  Pout  to  ko^  that  is  of  primary  importance  and 
not  the  individual  values. 


6  Discussion 


We  have  proposed  a  procedure  for  modelling  outliers  and  fertility  jumps  when  estimating 
treatment  effects  in  agricultural  field  trials.  As  a  first  stage,  we  have  specified  a  state-space 
modelling  framework  for  agricultural  trials  in  the  absence  of  such  anomalies  which  includes 
many  previously  proposed  models  as  special  ca.ses.  One  special  case  which  has  not  received 
much  previous  attention,  the  so-called  partial  second  differencing  model  when  A]  =  X^. 
seems  particularly  interesting  because  of  its  properties  and  its  success  in  examples.  The 
overall  procedure  seems  to  be  successful  in  accommodating  outliers  and  fertility  jumps  in  a 
routine  way,  and  it  also  allows  us  to  take  account  of  uncertainty  about  model  specification 
and  the  presence  of  anomalies  when  making  inference  about  the  treatment  effects.  The 
framework  allows  for  both  random  and  fixed  effects  modelling  in  a  natural  way,  and  permits 
the  incorporation  of  prior  information  very  easily. 


The  procedure  we  have  used  to  model  fertility  jumps  can  be  employed  for  any  stationary 
fertility  model  by  again  assuming  independence  between  fertilities  across  the  jump.  For 
the  first  difference,  or  random  walk,  fertility  model,  the  Ft  —  Ft-\  are  iid  N(0,<7^^jj),  and  we 
could  model  a  fertility  jump  at  plot  i  by  multiplying  the  variance  by  a  constant 
(=100  say).  While  this  has  the  desired  effect,  examples  suggest  that  after  accounting  for 
possible  fertility  jumps  a  smoother  fertility  model  is  likely  to  be  more  appropriate.  Modelling 
fertility  jumps  under  models  such  as  the  ARIMA  (1,1,0)  or  second  differencing  models  is  not 
so  straightforward. 

There  has  been  little  attempt  to  model  fertility  jumps  in  agricultural  field  trials,  although 
the  occurrence  of  level  shifts  in  general  time  series  is  often  recognised.  Outliers  are  usually- 
dealt  with  by  replacing  these  observations  with  missing  values.  The  decision  to  designate 
an  observation  as  an  outlier  may  be  made  by  a  subjective  examination  of  residuals,  or  by 
a  formal  statistical  test  (see  for  example  Kitagawa,  1979).  Note  that  this  is  approximately 
equivalent  to  conditioning  on  a  single  submodel,  which  we  have  shown  to  be  unsatisfactory 
in  some  circumstances. 

The  posterior  distribution  of  r  can  also  be  approximated  using  the  Gibbs  sampler  (Geman 
and  Geman,  19S4;  Gelfand  and  Smith,  1990).  Here  we  define  the  state  X  —  [F^ 
and  from  an  initial  value  of  A'  replace  A',  with  a  realisation  from  the  distribution  of  ( A,l.\j(j  ^ 
r),  K).  .All  of  these  distributions  are  trivial  except  perhaps  when  A,  is  an  element  of  0  (Taplin, 
1990).  After  enough  replacements,  X  will  be  approximately  a  realisation  from  the  posterior 
p{X\Y),  and  hence  a  stochastic  approximation  can  be  formed  by  such  repeated  sampling. 

This  procedure  is  efficient  at  determining  C,  as  it  visits  the  submodels  in  time  propor¬ 
tional  to  their  posterior  probability.  For  well  designed  agricultural  field  trials,  however,  it 
appears  relatively  ecisy  to  determine  C.  Furthermore,  the  Gibbs  sampler  is  not  as  efficient  at 
approximating  p(r|T,  Me)  as  the  approximation  p{T\y\  Mc,0)  used  in  equation  (6),  which  uses 
the  normality  assumptions  to  better  advantage.  Another  disadvantage  to  the  Gibbs  sampler 
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is  that  the  estimation  must  be  repeated  for  different  submodel  priors,  so  that  a  sensitix  ity 
analysis  suca  as  Figure  11  in  Section  5  is  not  as  immediately  available.  Also,  by  examining 
the  submodels  included  in  C,  the  data,  and  using  the  intuition  behind  what  the  submodels 
represent,  it  may  be  possible  to  determine  if  enough  submodels  ha\  e  been  included  in  C.  On 
the  other  hand,  the  partial  results  available  to  date  (Raftery  and  Lewis.  1991)  suggest  that 
the  Gibbs  sampler  would  have  to  be  run  for  a  large  number  of  iterations  (perhaps  on  the 
order  of  4,000-5,000)  to  obtain  accuracy  comparable  to  that  obtained  with  the  approxima¬ 
tions  we  have  used  here.  Thus,  the  Gibbs  sampler  may  well  be  computationally  inefficient 
compared  to  the  approximations  used  here. 

Based  on  our  examples  we  suggest  the  following  points: 

(1)  It  is  better  to  account  for  uncertainty  about  outliers  and  fertility  Jumps  than  to  con¬ 
dition  on  a  single  set  of  such  anomalies. 

(2)  Possible  outliers  and  fertility  jumps  can  have  a  greater  influence  on  our  assessment  of 
the  treatment  effects  than  the  choice  of  any  (reasonable)  fertility  model. 

(3)  Accounting  for  possible  outliers  and  fertility  jumps  can  reduce  the  disparity  between 
treatment  estimates  from  different  fertility  models. 

(4)  The  possibility  of  outliers  and  fertility  jumps  should  be  considered  xvhilt  deciding  on 
an  appropriate  fertility  model  (if  one  such  fertility  model  is  to  be  conditioned  on). 
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_ Table  1:  comparison  of  submodels  for  trial  ARC  8. _ 

model  submodel  with  C23=l  submodel  with  0^3=2 
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Table  2:  comparison  of  submodels  for  trial  ARC  6. 
model  submodel  with  C5=l  submodel  with  05=8 

_ loglike  cr^  Ai _ loglike  A] 

A2=0  0.6  0.36  0.94  2.8  0.50  0.97 

Ai=A2  0.0  0.30  0.67  4.4  0.36  0.72 


Table  3:  The  12  best  submodels  for  trial  ARC  1  and  their  approximate  posterior  probabilities 

submodel  posterior 
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Table  4:  Marginal  posterior  probabilities  of  outliers  axid  fertility  jumps  for  trial  ARC  1 
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Figure  2:  posterior  distribution  for  early  spraying  -  no  spraying  in  the  mildew  trial;  solid 
line  true  posterior,  dotted  line  approximation  (5). 


Figure  3:  Comparison  between  fixed  and  random  effects  models,  trial  ARC  2. 

Posterior  means  for  varieties  (left)  and  means  divided  by  standard  deviations  for  treatment 
differences  (right). 

The  vertical  and  horizontal  axes  indicate  random  and  fixed  effect  models  respectively. 
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Figure  4:  effect  on  posterior  probability  of  an  outlier  when  conditioning  on  the  MLE  for 
T  under  the  simple  model  with  F=0  and  <Tobs=l.  Only  the  observation  yl  in  question  is 
non-zero.  The  three  graphs  assume  (left  to  right)  3,5,  and  9  rephcates  respectively,  with  the 
true  posterior  (solid)  and  approximation  (dotted). 
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Figure  5:  effect  on  posterior  probability  of  an  outlier  when  conditioning  on  the  MLE  for 
(T^bs  under  the  simple  model  with  F=0.  Only  the  observation  yl  in  question  is  non-zero. 
The  graphs  (left  to  right)  correspond  to  trials  with  (treatment, replicates)  combinations  of 
(5,3),  (10,3)  and  (20,3)  respectively.  The  true  posterior  (and  approximation  (7))  is  solid  and 
approximation  (8)  is  dotted. 
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Figure  6;  trial  ARC  8,  treatment  adjusted  yield  graphed  with  symbol  treatment  number 
under  submodel  c<=l  VL  Lines  are  estimated  fertility  from  Ai=A2  model  (solid)  and  A2=0 
model  (dotted). 


yield 


Figure  7:  trial  ARC  8,  treatment  adjusted  yield  graphed  with  symbol  treatment  number 
under  submodel  with  C23=2,  Ct  =  \  V<^23.  Lines  are  estimated  fertility  from  Ai=A2  model 
(solid)  and  A2=0  model  (dotted). 
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Figure  8:  trial  ARC  6,  treatment  adjusted  yield  graphed  with  symbol  treatment  number 
under  submodel  C(=l  V<.  Line  is  estiinated  fertility  from  Ai=A2  model. 
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Figure  9:  trial  ARC  1,  treatment  adjusted  yield  graphed  with  symbol  treatment  number 
under  submodel  Ct=l  Vt.  Line  is  estimated  fertility  from  Ai=A2  model. 
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Figure  10:  trial  ARC  1,  treatment  adjusted  yield  graphed  with  symbol  treatment  number 
under  submodel  Cs—S,  c<=l  V<^5.  Line  is  estimated  fertility  from  Ai=A2  model. 
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Figure  11:  effect  on  posterior  for  a  treatment  contrast  of  Poui  and  Pjump-  Lines  represent 
posterior  mean  (solid)  and  2^%,97^%  (dotted),  and  ^%,99^%  (dashed)  quantiles.  Stars 
represent  values  when  pout  or  Pjump  —  100%.  Graphs  from  left  refer  to  trial  ARC  8  (A2=0), 
trial  ARC  6  (Ai=A2),  and  trial  ARC  6  (A2=0). 
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